eReefs Xarray interactiveΒΆ

import os
import numpy as np
import pandas as pd
import xarray as xr

import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
cartopy.config['data_dir'] = os.getenv('CARTOPY_DIR', cartopy.config.get('data_dir'))

import cmocean
import hvplot.xarray

import holoviews as hv
from holoviews import opts, dim

import geoviews as gv
import geoviews.feature as gf
from geoviews import tile_sources as gvts

from cartopy import crs

gv.extension('bokeh')

import warnings
warnings.filterwarnings("ignore", category=RuntimeWarning) 
year = 2018
base_url = "http://thredds.ereefs.aims.gov.au/thredds/dodsC/s3://aims-ereefs-public-prod/derived/ncaggregate/ereefs/gbr4_v2/daily-monthly/EREEFS_AIMS-CSIRO_gbr4_v2_hydro_daily-monthly-"
hydrofiles = [f"{base_url}{year}-{month:02}.nc" for month in range(3, 5)]
hydrofiles
['http://thredds.ereefs.aims.gov.au/thredds/dodsC/s3://aims-ereefs-public-prod/derived/ncaggregate/ereefs/gbr4_v2/daily-monthly/EREEFS_AIMS-CSIRO_gbr4_v2_hydro_daily-monthly-2018-03.nc',
 'http://thredds.ereefs.aims.gov.au/thredds/dodsC/s3://aims-ereefs-public-prod/derived/ncaggregate/ereefs/gbr4_v2/daily-monthly/EREEFS_AIMS-CSIRO_gbr4_v2_hydro_daily-monthly-2018-04.nc']
ds_hydro = xr.open_mfdataset(hydrofiles)
ds_hydro
<xarray.Dataset>
Dimensions:      (k: 17, latitude: 723, longitude: 491, time: 61)
Coordinates:
  * time         (time) datetime64[ns] 2018-02-28T14:00:00 ... 2018-04-29T14:...
    zc           (k) float64 dask.array<chunksize=(17,), meta=np.ndarray>
  * latitude     (latitude) float64 -28.7 -28.67 -28.64 ... -7.096 -7.066 -7.036
  * longitude    (longitude) float64 142.2 142.2 142.2 ... 156.8 156.8 156.9
Dimensions without coordinates: k
Data variables:
    mean_cur     (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 723, 491), meta=np.ndarray>
    salt         (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 723, 491), meta=np.ndarray>
    temp         (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 723, 491), meta=np.ndarray>
    u            (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 723, 491), meta=np.ndarray>
    v            (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 723, 491), meta=np.ndarray>
    mean_wspeed  (time, latitude, longitude) float32 dask.array<chunksize=(31, 723, 491), meta=np.ndarray>
    eta          (time, latitude, longitude) float32 dask.array<chunksize=(31, 723, 491), meta=np.ndarray>
    wspeed_u     (time, latitude, longitude) float32 dask.array<chunksize=(31, 723, 491), meta=np.ndarray>
    wspeed_v     (time, latitude, longitude) float32 dask.array<chunksize=(31, 723, 491), meta=np.ndarray>
Attributes: (12/21)
    Conventions:                     CF-1.0
    NCO:                             4.4.4
    Run_ID:                          2
    _CoordSysBuilder:                ucar.nc2.dataset.conv.CF1Convention
    aims_ncaggregate_buildDate:      2020-08-21T14:27:56+10:00
    aims_ncaggregate_datasetId:      products__ncaggregate__ereefs__gbr4_v2__...
    ...                              ...
    paramhead:                       GBR 4km resolution grid
    shoc_version:                    v1.1 rev(5620)
    technical_guide_link:            https://eatlas.org.au/pydio/public/aims-...
    technical_guide_publish_date:    2020-08-18
    title:                           eReefs AIMS-CSIRO GBR4 Hydrodynamic v2 d...
    DODS_EXTRA.Unlimited_Dimension:  time
reef_lat = -18.82
reef_lon = 147.64
min_lon = 146.5
min_lat = -20
max_lon = 148
max_lat = -17

lon_bnds = [min_lon, max_lon]
lat_bnds = [min_lat, max_lat]
ds_hydro_clip = ds_hydro.sel(latitude=slice(*lat_bnds), longitude=slice(*lon_bnds))
ds_hydro_clip.coords['k'] = ('zc',ds_hydro_clip.zc)
ds_hydro_clip = ds_hydro_clip.swap_dims({'zc':'k'})
ds_hydro_clip = ds_hydro_clip.drop(['zc'])
ds_hydro_clip
<xarray.Dataset>
Dimensions:      (k: 17, latitude: 100, longitude: 50, time: 61)
Coordinates:
  * time         (time) datetime64[ns] 2018-02-28T14:00:00 ... 2018-04-29T14:...
  * latitude     (latitude) float64 -20.0 -19.97 -19.94 ... -17.09 -17.06 -17.03
  * longitude    (longitude) float64 146.5 146.5 146.6 ... 147.9 148.0 148.0
  * k            (k) float64 -145.0 -120.0 -103.0 -88.0 ... -5.55 -3.0 -1.5 -0.5
Data variables:
    mean_cur     (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 100, 50), meta=np.ndarray>
    salt         (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 100, 50), meta=np.ndarray>
    temp         (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 100, 50), meta=np.ndarray>
    u            (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 100, 50), meta=np.ndarray>
    v            (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 100, 50), meta=np.ndarray>
    mean_wspeed  (time, latitude, longitude) float32 dask.array<chunksize=(31, 100, 50), meta=np.ndarray>
    eta          (time, latitude, longitude) float32 dask.array<chunksize=(31, 100, 50), meta=np.ndarray>
    wspeed_u     (time, latitude, longitude) float32 dask.array<chunksize=(31, 100, 50), meta=np.ndarray>
    wspeed_v     (time, latitude, longitude) float32 dask.array<chunksize=(31, 100, 50), meta=np.ndarray>
Attributes: (12/21)
    Conventions:                     CF-1.0
    NCO:                             4.4.4
    Run_ID:                          2
    _CoordSysBuilder:                ucar.nc2.dataset.conv.CF1Convention
    aims_ncaggregate_buildDate:      2020-08-21T14:27:56+10:00
    aims_ncaggregate_datasetId:      products__ncaggregate__ereefs__gbr4_v2__...
    ...                              ...
    paramhead:                       GBR 4km resolution grid
    shoc_version:                    v1.1 rev(5620)
    technical_guide_link:            https://eatlas.org.au/pydio/public/aims-...
    technical_guide_publish_date:    2020-08-18
    title:                           eReefs AIMS-CSIRO GBR4 Hydrodynamic v2 d...
    DODS_EXTRA.Unlimited_Dimension:  time
# Take the top 10 m
slice_ds = ds_hydro_clip.sel(k=slice(-10, 0)).drop(['u','v','eta','mean_wspeed','wspeed_u','wspeed_v'])
slice_ds
<xarray.Dataset>
Dimensions:    (k: 5, latitude: 100, longitude: 50, time: 61)
Coordinates:
  * time       (time) datetime64[ns] 2018-02-28T14:00:00 ... 2018-04-29T14:00:00
  * latitude   (latitude) float64 -20.0 -19.97 -19.94 ... -17.09 -17.06 -17.03
  * longitude  (longitude) float64 146.5 146.5 146.6 146.6 ... 147.9 148.0 148.0
  * k          (k) float64 -8.8 -5.55 -3.0 -1.5 -0.5
Data variables:
    mean_cur   (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 5, 100, 50), meta=np.ndarray>
    salt       (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 5, 100, 50), meta=np.ndarray>
    temp       (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 5, 100, 50), meta=np.ndarray>
Attributes: (12/21)
    Conventions:                     CF-1.0
    NCO:                             4.4.4
    Run_ID:                          2
    _CoordSysBuilder:                ucar.nc2.dataset.conv.CF1Convention
    aims_ncaggregate_buildDate:      2020-08-21T14:27:56+10:00
    aims_ncaggregate_datasetId:      products__ncaggregate__ereefs__gbr4_v2__...
    ...                              ...
    paramhead:                       GBR 4km resolution grid
    shoc_version:                    v1.1 rev(5620)
    technical_guide_link:            https://eatlas.org.au/pydio/public/aims-...
    technical_guide_publish_date:    2020-08-18
    title:                           eReefs AIMS-CSIRO GBR4 Hydrodynamic v2 d...
    DODS_EXTRA.Unlimited_Dimension:  time
vards = slice_ds.mean(dim='k').salt.load()
coastline = gf.coastline(line_width=3, line_color='k').opts(projection=ccrs.PlateCarree(), scale='10m')
land = gf.land.options(scale='10m', fill_color='lightgray')
var_mask = gv.Dataset(vards, crs=crs.PlateCarree())
minvar = vards.min().item()
maxvar = vards.max().item()
var0 = vards.fillna(0)

hv_ds = hv.Dataset(var0)

# Create stack of images grouped by time
im_mask = var_mask.to(gv.Image, ['longitude', 'latitude'], dynamic=True)
hv.output(widget_location='bottom')
label = slice_ds.salt.long_name+' '+slice_ds.salt.units
image = im_mask.opts(active_tools=['wheel_zoom', 'pan'], cmap=cmocean.cm.curl,
                     colorbar=True, width=450, height=400, clim=(34,36),
                     title=label) * coastline * land

image
/usr/share/miniconda/envs/envireef/lib/python3.8/site-packages/cartopy/io/__init__.py:260: DownloadWarning: Downloading: https://naciscdn.org/naturalearth/110m/physical/ne_110m_land.zip
  warnings.warn('Downloading: {}'.format(url), DownloadWarning)
opts.defaults(
    opts.GridSpace(shared_xaxis=True, shared_yaxis=True),
    opts.Image(cmap=cmocean.cm.curl, clim=(34,36)),
    opts.Labels(text_color='white', text_font_size='8pt', text_align='left', text_baseline='bottom'),
    opts.Path(color='white'),
    opts.Spread(width=600),
    opts.Overlay(show_legend=False))
im_hist = var_mask.to(gv.Image, ['longitude', 'latitude']).hist().opts(width=500, height=450, title=label)
im_hist
hv.save(im_hist, 'hist.html')
  0%|          | 0/61 [00:00<?, ?it/s]
  2%|▏         | 1/61 [00:00<00:18,  3.18it/s]
  3%|β–Ž         | 2/61 [00:00<00:18,  3.17it/s]
  5%|▍         | 3/61 [00:00<00:18,  3.18it/s]
  7%|β–‹         | 4/61 [00:01<00:17,  3.18it/s]
  8%|β–Š         | 5/61 [00:01<00:17,  3.18it/s]
 10%|β–‰         | 6/61 [00:01<00:17,  3.19it/s]
 11%|β–ˆβ–        | 7/61 [00:02<00:16,  3.18it/s]
 13%|β–ˆβ–Ž        | 8/61 [00:02<00:16,  3.17it/s]
 15%|β–ˆβ–        | 9/61 [00:02<00:16,  3.16it/s]
 16%|β–ˆβ–‹        | 10/61 [00:03<00:16,  3.16it/s]
 18%|β–ˆβ–Š        | 11/61 [00:03<00:15,  3.16it/s]
 20%|β–ˆβ–‰        | 12/61 [00:03<00:15,  3.16it/s]
 21%|β–ˆβ–ˆβ–       | 13/61 [00:04<00:15,  3.17it/s]
 23%|β–ˆβ–ˆβ–Ž       | 14/61 [00:04<00:14,  3.14it/s]
 25%|β–ˆβ–ˆβ–       | 15/61 [00:04<00:14,  3.09it/s]
 26%|β–ˆβ–ˆβ–Œ       | 16/61 [00:05<00:14,  3.06it/s]
 28%|β–ˆβ–ˆβ–Š       | 17/61 [00:05<00:14,  3.07it/s]
 30%|β–ˆβ–ˆβ–‰       | 18/61 [00:05<00:14,  3.06it/s]
 31%|β–ˆβ–ˆβ–ˆ       | 19/61 [00:06<00:13,  3.07it/s]
 33%|β–ˆβ–ˆβ–ˆβ–Ž      | 20/61 [00:06<00:13,  3.10it/s]
 34%|β–ˆβ–ˆβ–ˆβ–      | 21/61 [00:06<00:12,  3.12it/s]
 36%|β–ˆβ–ˆβ–ˆβ–Œ      | 22/61 [00:07<00:12,  3.14it/s]
 38%|β–ˆβ–ˆβ–ˆβ–Š      | 23/61 [00:07<00:12,  3.14it/s]
 39%|β–ˆβ–ˆβ–ˆβ–‰      | 24/61 [00:07<00:11,  3.15it/s]
 41%|β–ˆβ–ˆβ–ˆβ–ˆ      | 25/61 [00:07<00:11,  3.16it/s]
 43%|β–ˆβ–ˆβ–ˆβ–ˆβ–Ž     | 26/61 [00:08<00:11,  3.16it/s]
 44%|β–ˆβ–ˆβ–ˆβ–ˆβ–     | 27/61 [00:08<00:10,  3.16it/s]
 46%|β–ˆβ–ˆβ–ˆβ–ˆβ–Œ     | 28/61 [00:08<00:10,  3.16it/s]
 48%|β–ˆβ–ˆβ–ˆβ–ˆβ–Š     | 29/61 [00:09<00:10,  3.15it/s]
 49%|β–ˆβ–ˆβ–ˆβ–ˆβ–‰     | 30/61 [00:09<00:09,  3.15it/s]
 51%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆ     | 31/61 [00:09<00:09,  3.16it/s]
 52%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–    | 32/61 [00:10<00:09,  3.16it/s]
 54%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–    | 33/61 [00:10<00:08,  3.13it/s]
 56%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–Œ    | 34/61 [00:10<00:08,  3.12it/s]
 57%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–‹    | 35/61 [00:11<00:08,  3.13it/s]
 59%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–‰    | 36/61 [00:11<00:07,  3.14it/s]
 61%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆ    | 37/61 [00:11<00:07,  3.14it/s]
 62%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–   | 38/61 [00:12<00:07,  3.15it/s]
 64%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–   | 39/61 [00:12<00:06,  3.16it/s]
 66%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–Œ   | 40/61 [00:12<00:06,  3.16it/s]
 67%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–‹   | 41/61 [00:13<00:06,  3.16it/s]
 69%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–‰   | 42/61 [00:13<00:06,  3.17it/s]
 70%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆ   | 43/61 [00:13<00:05,  3.18it/s]
 72%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–  | 44/61 [00:13<00:05,  3.18it/s]
 74%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–  | 45/61 [00:14<00:05,  3.18it/s]
 75%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–Œ  | 46/61 [00:14<00:04,  3.18it/s]
 77%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–‹  | 47/61 [00:14<00:04,  3.17it/s]
 79%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–Š  | 48/61 [00:15<00:04,  3.16it/s]
 80%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆ  | 49/61 [00:15<00:03,  3.17it/s]
 82%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ– | 50/61 [00:15<00:03,  3.16it/s]
 84%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–Ž | 51/61 [00:16<00:03,  3.16it/s]
 85%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–Œ | 52/61 [00:16<00:02,  3.17it/s]
 87%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–‹ | 53/61 [00:16<00:02,  3.16it/s]
 89%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–Š | 54/61 [00:17<00:02,  3.17it/s]
 90%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆ | 55/61 [00:17<00:01,  3.17it/s]
 92%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–| 56/61 [00:17<00:01,  3.17it/s]
 93%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–Ž| 57/61 [00:18<00:01,  3.17it/s]
 95%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–Œ| 58/61 [00:18<00:00,  3.11it/s]
 97%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–‹| 59/61 [00:18<00:00,  3.07it/s]
 98%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–Š| 60/61 [00:19<00:00,  3.07it/s]
100%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆ| 61/61 [00:19<00:00,  3.06it/s]
                                               

slice_ds.coords['y'] = ('latitude',slice_ds.latitude)
slice_ds.coords['x'] = ('longitude',slice_ds.longitude)
slice_ds2 = slice_ds.swap_dims({'latitude':'y'})
slice_ds2 = slice_ds2.swap_dims({'longitude':'x'})
slice_ds2 = slice_ds2.drop(['latitude'])
slice_ds2 = slice_ds2.drop(['longitude'])
slice_ds2 =  slice_ds2.assign_coords( x=(np.arange(slice_ds.x.shape[0])),
                                 y=(np.arange(slice_ds.y.shape[0])),
                                 ) 
slice_ds2
<xarray.Dataset>
Dimensions:   (k: 5, time: 61, x: 50, y: 100)
Coordinates:
  * time      (time) datetime64[ns] 2018-02-28T14:00:00 ... 2018-04-29T14:00:00
  * k         (k) float64 -8.8 -5.55 -3.0 -1.5 -0.5
  * y         (y) int64 0 1 2 3 4 5 6 7 8 9 10 ... 90 91 92 93 94 95 96 97 98 99
  * x         (x) int64 0 1 2 3 4 5 6 7 8 9 10 ... 40 41 42 43 44 45 46 47 48 49
Data variables:
    mean_cur  (time, k, y, x) float32 dask.array<chunksize=(31, 5, 100, 50), meta=np.ndarray>
    salt      (time, k, y, x) float32 dask.array<chunksize=(31, 5, 100, 50), meta=np.ndarray>
    temp      (time, k, y, x) float32 dask.array<chunksize=(31, 5, 100, 50), meta=np.ndarray>
Attributes: (12/21)
    Conventions:                     CF-1.0
    NCO:                             4.4.4
    Run_ID:                          2
    _CoordSysBuilder:                ucar.nc2.dataset.conv.CF1Convention
    aims_ncaggregate_buildDate:      2020-08-21T14:27:56+10:00
    aims_ncaggregate_datasetId:      products__ncaggregate__ereefs__gbr4_v2__...
    ...                              ...
    paramhead:                       GBR 4km resolution grid
    shoc_version:                    v1.1 rev(5620)
    technical_guide_link:            https://eatlas.org.au/pydio/public/aims-...
    technical_guide_publish_date:    2020-08-18
    title:                           eReefs AIMS-CSIRO GBR4 Hydrodynamic v2 d...
    DODS_EXTRA.Unlimited_Dimension:  time
saltXY = slice_ds2.mean(dim='k').salt.load()
hv_ds = hv.Dataset(saltXY)

# Create stack of images grouped by time
im = hv_ds.to(hv.Image, ['x','y'], dynamic=True).opts(active_tools=['wheel_zoom', 'pan'], cmap=cmocean.cm.curl,
                     colorbar=True, width=450, height=400, clim=(34,36))

polys = hv.Polygons([])

box_stream = hv.streams.BoxEdit(source=polys)

# Declare an empty DataFrame to declare the types
empty = pd.DataFrame({'time': np.array([], dtype='datetime64[ns]'), 'salt': []})

def roi_curves(data):
    if not data or not any(len(d) for d in data.values()):
        return hv.NdOverlay({0: hv.Curve(empty, 'time', 'salt')})

    curves = {}
    data = zip(data['x0'], data['x1'], data['y0'], data['y1'])
    for i, (x0, x1, y0, y1) in enumerate(data):
        selection = hv_ds.select(x=(x0, x1), y=(y0, y1))
        curves[i] = hv.Curve(selection.aggregate('time', np.mean))
    return hv.NdOverlay(curves)

# Generate VLines by getting time value from the image frames
def vline(frame):
    return hv.VLine(frame.data.time.values)
vlines = im.apply(vline)

dmap = hv.DynamicMap(roi_curves, streams=[box_stream])

(im * polys + dmap * vlines ).opts(
    opts.Curve(width=400, framewise=True), 
    opts.Polygons(fill_alpha=0.2, line_color='white'), 
    opts.VLine(color='black')).opts(title=label)